■■nfi 


um ' 


CENTRAL  CIRCULATION  BOOKSTACKS 

The  person  charging  this  material  is  re- 
sponsible for  its  renewal  or  its  return  to 
the  library  from  which  it  was  borrowed 
on  or  before  the  Latest  Date  stamped 
below.  You  may  be  charged  a  minimum 
fee  of  $75.00  for  each  lost  book. 

Theft,  mutilation,  and  underlining  of  books  are  reasons 
for  disciplinary  action  and  may  result  in  dismissal  from 
the  University. 

TO  RENEW  CALL  TELEPHONE  CENTER,  333-8400 
UNIVERSITY    OF    ILLINOIS    LIBRARY   AT    URBANA-CHAMPAIGN 


JUN  1  8  1936 
SEP  23i997 
JUN  0  3  1997 


When  renewing  by  phone,  write  new  due  date  below 
previous  due  date.  L162 


Digitized  by  the  Internet  Archive 
in  2013 


http://archive.org/details/numericalintegra894fatu 


7  £dg 

'U^U&CDCS-R-77-894 


COO-2383-0044 

UILU-ENG  77  1751 


^7 


NUMERICAL  INTEGRATORS  FOR  STIFF 
AND  HIGHLY  OSCILLATORY  DIFFERENTIAL  EQUATIONS 


by 


Simeon  Ola  Fatunla 


October  1977 


UIUCDCS-R-77-894 


NUMERICAL  INTEGRATORS  FOR  STIFF 
AND  HIGHLY  OSCILLATORY  DIFFERENTIAL  EQUATIONS 

by 

Simeon  Ola  Fatunla* 


October  1977 


DEPARTMENT  OF  COMPUTER  SCIENCE 

UNIVERSITY  OF  ILLINOIS  AT  URB ANA- CHAMPAIGN 

URBANA,  ILLINOIS  61801 


*Author  on  leave  from  the  Department  of  Mathematics,  University  of  Benin, 
Benin-City,  Nigeria.   Work  supported  in  part  by  the  U.  S.  Energy  Research 
and  Development  Administration  under  contract  US  ERDA/EY-76-S-02-2383. 


ACKNOWLEDGMENT 


The  author  is  indebted  to  Professor  C.  W.  Gear  for  his 
invaluable  suggestions,  criticism  and  advice.   He  has  gone  out  of  his 
way  to  ensure  that  I  was  comfortable  during  my  short  visit  at  the 
Department  of  Computer  Science. 


-1- 


ABSTRACT 


Some  L-stable  fourth  order  explicit  one  step  numerical  integration 
formulas  which  require  no  matrix  inversion  are  proposed  to  cope  effectively 
with  systems  of  ordinary  differential  equations  with  large  Lipschitz 
constants  (including  those  having  highly  oscillatory  solutions).   The 
implicit  integration  procedure  proposed  in  Fatunla  [10]  is  further  developed 
to  handle  a  larger  class  of  stiff  systems  as  well  as  those  with  highly 
oscillatory  solutions.   The  same  pair  of  nonlinear  equations  as  in  [10]  is 
solved  for  the  stiffness/oscillatory  parameters.   However,  the  nonlinear 
systems  are  transformed  into  linear  forms  and  an  efficient  computational 
procedure  is  developed  to  obtain  these  parameters.   The  new  schemes  compare 
favorably  with  the  backward  differentiation  formula  (DIFSUB)  of  Gear  [12, 
13]  and  the  blended  linear  multistep  methods  of  Skeel  and  Kong  [22],  and 
the  symmetric  multistep  methods  of  Lambert  and  Watson  [16]. 
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1.   INTRODUCTION 

The  development  of  numerical  integration  formulas  for  stiff  as 
well  as  highly  oscillatory  systems  of  differential  equations  has  attracted 
considerable  attention  in  the  past  decade.   The  reason  for  this  cannot  be 
farfetched,  realizing  that  the  mathematical  models  of  physical  situations 
in  kinetic  chemical  reactions,  process  control  and  electrical  circuit 
theory  often  generate  systems  of  ordinary  differential  equations  whose 
Jacobians  have  at  least  one  eigenvalue  with  a  very  negative  real  part  or 
very  large  imaginary  part.   Both  situations  are  respectively  described 
as  stiff  and  highly  oscillatory. 

Consider  the  following  model  problems 
(1.1)      y'  =  A(y-c(x))  +  c'(x), 

y(a)  =  yQ, 
where  y(x)  £  R  ,  X   a  complex  constant  with  re  X   <<  0  and  c(x)  has  the 
property  |c'(x)|  ~   0  in  the  finite  interval  a  <  x  <  b; 


(1.2)  !-e      J, 


z'  =        z»    z(a)  =  zQ, 

'  -co  -e 


2 
with  z(x)  E  R  and  co  >>  0,  e  a  positive  constant  such  that  e    ~    0. 

Problem  (1.1)  has  theoretical  solution 
(1.3)     y(x)  =  c(x)  +  yQeAx 

whose  component  c(x)  is  slowly  varying  in  the  specified  interval  while  the 
second  component  decays  rapidly  in  the  transient  phase. 

The  analytic  solution  to  problem  (1.2)  is  given  by 


(1.4)  /  a  cos  cox  +  6  sin  cox'. 

z(x)  =  e  CX 

1+8  cos  cox  -  a  sin  cox, 

where  a  and  6  are  the  arbitrary  constants  of  integration. 
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The  transitory  phase  for  problem  (1.1)  is  of  the  order  of  -1/A 
while  that  of  problem  (1.2)  is  the  entire  interval  of  integration  with 
o)/2tt  complete  oscillations  per  unit  interval. 

Almost  invariably,  most  conventional  numerical  integration 
solvers  cannot  effectively  cope  with  problems  (1.1)  and  (1.2)  as  they  lack 
adequate  stability  characteristics.   Any  attempt  to  impose  the  stability 
properties  will  in  effect  constrain  the  integration  mesh-size  to  be 
intolerably  small.   This  may  ultimately  have  adverse  effects  on  the 
accuracy  due  to  propagation  of  roundoff  errors.   Besides,  the  computing 
time  and  cost  may  be  too  excessive. 

Existing  algorithms  developed  for  problems  of  the  type  (1.1) 
can  be  classified  into  the  following  categories: 

(i)   generalized  Runge-Kutta  schemes 
(ii)   implicit  Runge-Kutta  schemes 
(iii)   trapezoidal  rule  with  extrapolation 
(iv)   multiderivative  multistep  formulas 

(v)   backward  differentiation  multistep  formulas. 

As  of  now  the  most  widely  used  numerical  integration  code  for 
stiff  systems  is  the  class  (v)  schemes  particularly  Gear's  DIFSUB  [12,  13]. 
The  method  is  efficient  and  reliable  provided  the  eigenvalues  of  the 
Jacobian  are  not  close  to  the  imaginary  axis  where  the  higher  order  schemes 
exhibit  poor  stability  properties  (as  evidenced  from  example  2).   Dahlquist 
[5]  established  that  neither  an  explicit  linear  multistep  scheme  of  any 
order  nor  an  implicit  multistep  method  of  order  greater  than  two  can  be 

A-stable.   He  proved  that  the  trapezoidal  rule  (which  of  course  is  not 

1   3  (3) 

L-stable  as  y  /y   ,  -*■   1  as  Ah  -*  -°°)  has  the  smallest  error  of  (±)to  h  y.  ,. 
Jt      t-1  1^     V.x; 

Other  schemes  which  behave  better  than  DIFSUB  when  the  eigenvalues  are  close 
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to  the  imaginary  axis  includes  the  second  derivative  multistep  method 
Enright  [6]. 

As  regards  problems  of  the  type  (1.2),  earlier  efforts  include  Gautschi'tl 
[11]  nonlinear  multistep  schemes  which  produce  exact  solution  to  algebraic 
(or  trigonometric)  polynomials  up  to  certain  degrees.   The  main  drawback 
of  this  scheme  is  that  it  requires  an  a  priori  knowledge  of  the  period  of 
the  systems  under  consideration.   Other  numerical  integration  solvers  to 
oscillatory  systems  include  Amdursky  and  Ziv  [1,  2,  3],  Snider  and 
Flemming  [23],  Miranker,  et  al.  [18],  Miranker  and  Wahba  [19],  Miranker 
and  Veldhuizen  [21],  and  Fatunla  [8,  9].   Unfortunately,  none  of  these 
existing  routines  has  been  properly  (adequately)  put  to  test  for  various 
kinds  of  oscillatory  problems. 

In  this  paper,  we  propose  some  two-point  numerical  integration 
formulas  which  effectively  cope  with  systems  of  ODEs  whose  characteristics 
are  identical  to  problems  (1.1)  and  (1.2). 

We  shall  consider  initial  value  problems 

(1.5)  y'  =  f(x,y),  y(0)  =  yQ 

with  y(x)  £  R  in  the  finite  interval  S  =  [0,xf]  C  R  where  x  =  Nh  for 
some  positive  integer  N  >  0.   It  is  assumed  that  y(x)  is  sufficiently 
dif ferentiable.   We  adopt  the  vector  notation:   y  =  (  y'   y, . . . ,   y)  , 
f=(f,   f,...,   f).   The  numerical  estimates  y   to  the  theoretical 
solution  y(x  )  at  the  points  x  =  th,  t  =  0(1)N  are  to  be  generated. 

On  every  subinterval  [x  ,  x  +h],  the  theoretical  solution  y(x) 
is  approximated  by  either  the  interpolating  function 

fix         -fi„x 

(1.6)  F(x)  =  (I-e    )A  -  (I-e  l    )B  +  C, 

A,  B,  nnd  C  being  m-tuples  with  real  entries,  I  is  the  identity  matrix 
whilst  ru     and  U     are  diagonal  (stiffness/oscillatory)  matrices  or 
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a  x       n*x 

(1.7)  F(x)  =  (I-e   )R  +  (I-e   )R*  +  S 

where  R,  S  are  m-tuples  with  complex  entries  and  (*)  denotes  complex 
conjugate. 

The  choice  of  interpolation  formula  is  determined  by  equation 
(3.12). 

The  following  definitions  are  worthwhile: 
Definition  1.   A  one-step  numerical  integration  scheme  is  considered  L-stable 
if  apart  from  being  A-stable,  when  it  is  applied  to  the  scalar  initial 
value  problem 

(1.8)  y'  =  Ay,  y(0)  =  n 

(X  being  a  complex  constant  with  negative  real  part) ,  the  resultant 
numerical  solution  is  given  by 

(1.9)  yfc+1  =  u(Xh)yt 

with  the  characteristic  equation  y(Xh)  having  the  property: 

(1.10)  lim   |y(Xh)|  =  0. 
re(Xh)-x» 

Definition  2.   A  numerical  integration  scheme  is  said  to  be  exponentially 

fitted  at  a  complex  value  X  =  X   if  when  it  is  applied  to  the  initial 

value  problem  (1.8)  with  exact  initial  condition,  the  characteristic 

equation  y (Xh)  satisfies  the  relation 

X  h 

(1.11)  y(XQh)  =  e  U  . 

Liniger  and  Willoughby  [18]  and  Jackson  and  Kenue  [14]  have 
respectively  discussed  A-stable  one-  and  two-step  numerical  integration 
methods  which  are  exponentially  fitted  at  infinity.   Both  schemes  ensure 
exponential  fitting  by  a  suitable  choice  of  a  free  parameter.   This 
approach  was  further  extended  to  construct  a  stiffly  stable  k-step  method 
of  order  k+2  in  Enright  [6]. 
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We  shall  construct  from  both  equations  (1.6)  and  (1.7)  explicit 
one-step  numerical  integration  formulas  of  fixed  order  four  which  possess 
adequate  stability  and  convergent  characteristics  to  cope  with  both  stiff 
and  highly  oscillatory  systems  of  ordinary  differential  equations.   For 
linear  systems  the  interpolating  functions  are  global  in  the  sense  that  the 
stiffness/oscillatory  matrices  have  constant  entries  which  are  determined 
by  solving  a  set  of  linear  equations  at  the  first  step  of  the  integration 
procedure.   The  implicit  scheme  proposed  in  Fatunla  [10]  is  further 
developed  using  the  interpolating  function  (1.7).   The  need  to  solve 
nonlinear  equations  for  the  stiffness/oscillatory  parameters  is  eliminated. 
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2.   DEVELOPMENT  OF  THE  INTEGRATION  FORMULAS 


Let  y  . .  denote  the  numerical  estimate  of  the  theoretical  solution 

t+j 

y(x  ,.)  and  x  =  x  ,  ,  and  adopt  the  notation  f  ,.  ■  f(x_..,  y_. .)  and  set 
7   t+j  t+j        r  t+j      t+j  Jt+3 

(2.1)  yt+.  ■  F(V]),    3=0,1 
and 

(2.2)  f^k)  =  i™y  k  =  0,  1. 

The  imposition  of.  these  constraints  on  the  interpolating  function 
(1.6)  gives  the  integration  formula 

(2.3)  yt+1  -  yt  +  «t  +  sf  <»  , 
where 

(2.4)  r  =  Oy-^a)  , 
and 

(2.5)  s  =  y  +o  ; 

Y  and  a  are  diagonal  matrices  with  entries 

xft  h 

(2.6)  iY=| r^ ,  i  =  l(l)m, 

1J21(1fi1+1fl2) 

and 

(2.7)  ^  =  -^ : t-^-  ,  i  =  l(l)m. 

1fi2(1fi1+1fi2) 

In  the  event  that  a  component  of  the  stiffness/oscillatory  matrix 
(ft  say)  does  vanish,  by  L'hopital  rule,  the  corresponding  component  of 

Y  in  (2.4)  is 

(2.8)  3    -h_ 

Jft 
"2 

and  the  resultant  integration  formula  is  given  by 
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(2.9)  jy    =  jy  +  h j  f  +  (-7^  +  ja)jf^1)  . 

We  now  consider  the  case  of  the  complex  interpolating  function 
(1.7).   Let  the  parameters  be  expressed  as 

(2.10)  tt     =  A  +  iu 


and 


* 

0,     =  A  -  iu 


By  imposing  constraints  (2.1)  and  (2.2)  on  (1.7),  we  obtained 
the  integration  formula  (2.3)  with  r  and  s  given  by 

,~   , ,  s      ,,   >.    e  {(A  -u  )sin(hu)-2Au  cos(hu)}  +  2Au 
\Z.A.L)  r^A,u;  -  *      r  , 

u(A  +u  ) 

\  v» 
/0  10N      ,,   v    e   {A  sin(hu)-u  cos(hu)}  +  u 

u(AZ+u  ) 
Near  the  origin,  equations  (2.11)  and  (2.12)  reduce  to 

(2.13)  r(A,u)  =  h 
and 

(2.14)  s(A,u)  =  ^-  , 

and  thus  the  integration  formula  (2.3)  reduces  to  the  second  order 
Taylor  series. 
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3.   EVALUATION  OF  STIFFNESS  PARAMETERS 


By  using  the  Taylors  expansion  of  y  ..  =  y(x  +h)  about  x  =  x 

«1h   -Oh  t 

and  the  Maclaurin's  series  of  e    ,  e     in  the  integration  formula  (2.3) 

0   1       2 
it  is  observed  that  the  coefficients  of  h  ,  h  ,  and  h  vanish  identically. 

With  the  view  to  obtain  numerical  estimates  for  the  stiffness/oscillatory 

3      4 
parameters,  we  simply  allow  the  coefficients  of  h  and  h  to  vanish  thus 

yielding  the  following  pair  of  nonlinear  equations 

(3.1)  (n2-fi1)f^1)  -  n1n2ft  -  -f£2), 

and 

(3.2)  -(flJ-fl^+G*)^  +  ^2(fi2-fi1)f^1)  =  -f^3)  . 

By  adopting  (3.1)  as  a  definition  of  £Lft  f  ,  equation  (3.2) 
becomes 

(3.3)  (fi2-fi1)f^2)  -  fi1"2^1)  =  ~f^3)  • 

If  the  set  of  equations  (3.1)  and  (3.3)  were  to  be  meaningful, 

it  is  desirable  that 

if(2)   i,(l) 
r      r 
t      t 

(3.4)  det  i  #  0,    i  =  l(l)m. 


\  I 


Let 


(3.5)  XD  =  X^2  -  lQ  ,    i  =  l(l)m; 
and 

(3.6)  iE  =  1fi11fi2  ,    i  =  l(l)m 

Equations  (3.1)  and  (3.3)  can  now  be  expressed  as  m  pairs  of 
linear  equations 


(3.7) 
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^Vf(2)  -  Vf(1)  =-V3)  , 
t        t        t   ' 

for  i  =  l(l)m. 


lDif(D  .  Vft   =  -if<2)  , 

These  m  pairs  of  equations  can  be  readily  solved  for  D  and 

E  by  the  Cramer's  rule  to  give 

if  if(3)  _  if(l)if(2) 

(3'8)     lD  ■  if(l)if(l)  .  I  if(2)   -    1=1(1>m 
t    t       t   t 

and 

if(l)if(3)  _  if(2)if(2) 

<3'9>    ±E  ■  if(i)if(i)  _  iffcif(2)  >        f"»-- 

t    t       t   t 
The  numerical  values  obtained  for  D  and  E  can  be  substituted 
in  equations  (3.5)  and  (3.6)  to  generate  the  oscillatory/stiffness 
parameters  as 


(3.10)    V  =  \[-\   ♦  /(V+^E)], 
and 

(3.ii)    xsi    =  xn  +  XD. 

The  complex  interpolation  formula  (1.7)  is  adopted  if  in 


(3.10),  the  following  relationship  holds 
(3.12)    V  <  4XE 
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4.   STABILITY  CONSIDERATIONS 

We  now  apply  the  integration  formula  (2.9)  to  the  scalar  test 
equation 

(4.1)  y'  =  Ay  +  g 

where  X  is  a  complex  constant  with  negative  real  part  and  g  is  a 
real  constant. 

The  numerical  solution  is  given  as 

(4.2)  yt+1  =  pU,fi2,h)yt  +  gh 

where  the  characteristic  equation  p(A,fl„,h)  is  given  by 

2      "Sh     2 

(4.3)  p(X,fl2,h)  =  1  +  Ah  +  X  (SI  h+e  -l)/«2 

Equation  (3.11)  gives  the  following  relationship: 

(4.4)  tt2   =  XD  =  D 

If  we  use  equation  (4.1)  in  (3.8),  we  readily  obtain 

(4.5)  D  =  -X  =  ft2 

This  in  equation  (4.3)  gives 

(4.6)  P(X,-X,h)  =  eAh  . 

We  further  consider  the  application  of  the  integration  formula 
(2.3)  to  test  problem  (4.1)  for  the  case  when  the  stiffness/oscillatory 
parameter  is  imaginary.   Thus,  r(X,u)  and  s(X,u)  used  in  the  two 
integration  formulas  (2.3)  to  the  test  problem  (4.1)  as  specified  by 

,,  ,v      fy      s        sin(hu) 

(4.7)  r(X,u)  =  J — '-   , 

and 

/,    Q\  /-,   n    l-cos(hu) 

(4.8)  s(X,u)  =  2 . 

u 


-12- 


Let 

(4.9)  X  =  iz,    i2  =  -1. 

The  resultant  integration  formula  is  given  by 

(4.10)  yt+1  =  q(X,u,h)yt  +  gh, 

where  the  characteristic  equation  can  be  obtained  as 

(4.11)  q(X,u,h)  =  1  +  (l-cos(hu))X2  +  sin(hu)X   < 

u 

The  application  of  equation  (4.1)  and  (4.9)  in  (3.8)  yields 

(4.12)  D  -  -iz  =  -X. 

■k 

=  n  . 

This  in  (4.11)  gives 

(4.13)  q(u,u,h)  =  1  -  (l-cos(hu))+i  sin(hu) 

ihu 
=  e 

Xh 
=  e 

Equations  (4.6)  and  (4.13)  together  with  definitions  1  and  2 

thus  establish  the  L-stability  and  exponential  fitting  of  the  proposed 

integration  formulas. 
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5.   LOCAL  TRUNCATION  ERROR 

We  now  associate  with  the  integration  formula  (2.3)  the  operator 
V[y(x),h]  specified  as 

(5.1)  V(y(x),h)  =  y(x+h)  -  y(x)  +  (Sl^-Qf  )f  (x,y)  -  (Y-to  )f (1)  (x,y) 

for  an  arbitrary  function  y(x)  G  c  (S).   The  local  truncation  error  T 
at  x  =  x  i  is  hence  given  as  V[y(x  ),h]  where  y(x  )  denotes  the  solution 
to  the  initial  value  problem  (1.5).   By  using  the  Taylors  expansion  of 
V[y(x),h]  about  x  =  x  with  the  localizing  assumption  that  there  is  no 
previous  error  (i.e.,  y  =  y(x  )),  the  truncation  error  for  the 

integration  formula  (2.3)  with  constraints  (3.1)  and  (3.3)  can  be  derived  at 

(5.2)  Tt+1  =  yj-{(fi1+^2)"1[(ni+fi2)f^4)  -  ^(ft(1)+fi2ft) 

+  ^2(f^1)-fi1ft)}+  0(h6). 

The  corresponding  truncation  formula  for  the  integration  formula 
(2.9)  is  given  by 

(5.3)  Tt+1=^  (f^-^a))  +0(h6}) 

while  the  truncation  error  when  r(A,u)  and  s(X,u)  are  specified  by  equations 
(2.11)  and  (2.12)  is 

(5.4)  Tt+1  =  3y  [f^4)  +  (A2+u2)(3A2-u2)ft  +  4X  (u2-A2)f^1)]  +  0(h6). 

From  equations  (5.3)  and  (5.4)  we  deduce  that  all  the  proposed 
explicit  integration  formulas  are  of  fixed  order  four. 
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6.   EXTENSION  OF  NONLINEAR  SCHEME 


We  recall  the  implicit  numerical  integration  formula  proposed 
in  Fatunla  [10] : 

(f,     -n  ,JS+1]     ■-    „tS]  TT         T[s]  n-1      r„[s]         f[s]  ,  ft        -         , 

(6.D    yt+1   -  yt+1  -  [i-  Jt+1l   [yt+1-  ft+rGt]'   s  "  °'  l'  2'" 

[si 
where  J.+1  denotes  the  Jacobian  specified  as 

(6  2)     J[s]  =  —  (x    v[s])   f[s]  =  f(x     v[s]) 

and 

(6.3)  Gt  =  yt  -  (Y-HT)ft. 

The  components  of  9,  y,  a  are  obtained  as 

ln  h  .  -ln.Ji 

Xfi    (1-e  )   +     fl..  (1-e        2   ) 

(6.4)  S-  — : ~ ,  i  =  l(l)m, 

fi     ft    (e  -e  ) 

-*Q  h  ln  h 

(6.5)  \=— (1"e   ,      }      ,  i  =   l(l)m, 

-   ft„h        flnh 
^1(e  -e  ) 

and 

Xft  h  _1^2h 

(6.6)  ^    =  — P^ : *-   ,  i  =   l(l)m. 

-Vh     Vh 
ft~(e  -e  ) 

The  corresponding  truncation  error  for  the  integration  formula 
(6.1)  was  obtained  as 

(6.7)  Tt+1  =  j^Q   {9ft4)  +  t(fi2"fil)ft3>  +  5("1-^2+fi2)ft2) 

+   [5n^«2  "  9n  n  (n^-n  n  -k^]^  +  °(h7)-  ■ 
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In  the  case  when  the  oscillatory/stiffness  parameters  have  complex 
values,  the  components  of  0  and  (y+<$)  in  equation  (6.1)  are  replaced  by 

re  o\  iQ  _  e   [u  cos(hu)  -  X  sin(hu)]  +  u       _  ,>.» 

(    }  ^7X~2~~7~r~r — ■       "  () 

e   (X  +u  )  sxn(hu) 

,,   nx     i,  .  ~N    ue   -  [X  sin(hu)  +  u  cos(hu)] 

(6.9)  (Y+6)  =  2 — 2 

(X  +u  )  sin(hu) 

with  the  truncation  error  (6.7)  replaced  by 

(6.10)  Tt+1  =  -jjQ   {9f^4)  -  lOXf <3)  +  5(3X2-u2)f^2) 

+  36(X2-u2)f^1)  +  4(X2+u2)(8X2-u2)ft}  +  0(h?)  . 
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NUMERICAL  EXAMPLES 


Example  1 


We  first  consider  the  linear  system: 


r 


-.1  -49.9 


(7.1) 


0 

y'  =  I  0   -50      0 

L  0    70    -120J 


y(0)  =!  l 

in  the  interval  0  £  x  _<  15.   The  eigenvalues  of  the  Jacobian  of  (7.1)  are 
^   =  -0.1,  ^?   =  -50,  ^„  =  -120.   The  theoretical  solution  is  given  as 

1  ,   x    -O.lx  ,   -50x 
y(x)  =  e      +  e 

2  ,  ,     -50x 
y(x)  =  e 

3  ,    s  -50x  ,   -120x 
y(x)  =  e     +  e 

The  new  explicit  scheme  performs  better  than  the  DIFSUB  [13]  and 

the  blended  DIFSUB  [22].   Details  of  the  numerical  results  are  given  in 

Table  7.1. 


TABLE  7.1.   Numerical  Results  for  Example  1. 

.,  .  ,    Max   _      Function   Back     LU    Accurate 
Method   _  ,     Steps    _  -.      _  -      _        _. :  .. 
Order  Calls    Solves  Decomp   Digits 


DIFSUB 

BLENDED 
DIFSUB 

EXPLICIT 

SCHEME 


6    353     1024      969     18 
11    234      595     1056     22 


75 


75 


9.3 

10.4 

12.5 
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Example  2 


We  further  consider  the  linear  problem: 


(7.2) 


~-10 

100 

0 

0 

0 

on 

-100 

-10 

0 

0 

0 

0 

0 

0 

-4 

0 

0 

0 

y.  = 

0 

0 

0 

-1 

0 

0 

0 

0 

0 

0 

-.5 

0 

0 

0 

0 

0 

-.1- 

y.(0)  =1,    i  =  1(1)6 


in  the  range  0  <_  t  <   20.   The  eigenvalues  of  the  Jacobian  are  A   _  =  -10  ±  lOOi, 

■•-»  * 

A-  =  -4,  A,  =  -1,  A,.  =  -.5,  and  A,  =  -.1.   This  problem  is  particularly 
troublesome  for  both  DIFSUB  and  the  blended  DIBSUB  as  can  be  seen  from  Table  7.2, 


TABLE  7.2.   Numerical  Results  for  Example  2. 


G 

Max           Func    Back 
Order          Eval   Solves 

LU 

Accurate 

Decomp 

Digits 

DIFSUB 

io"2 

(4)   (1001)   (3002)   (2959) 

(7) 

(1.1) 

"too  much  work" ;  integration  abandone 

.d  at  x  = 

8.3 

blended  DIFSUB 

io"10 

(12)   (1001)   (2917)   (5580) 

(21) 

(10.3) 

"too  much  work";  integration  abandone 

d  at  t  = 

5.1 

Explicit  Scheme 

unspec- 
ified 

4     200     200 

- 

14.2 
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Example  3 

We  now  consider  the  initial  value  problem  of  Liniger  and 
Willoughby  [18]: 


(7.3a) 


y'  = 


'-2000  iooo\    /i        jo) 
y  +  |   ,  y(0)  = 

\   0  ;  0 


\°/ 


in  the  interval  0  <  x  <  5. 

This  problem  is  contained  in  the  test  batch  recommended  by 
Bjurel,  et  al.  [4],  and  was  also  considered  in  Lambert  [15].   In  addition, 
it  is  solved  with  the  blended  DIFSUB  [22]. 

The  eigenvalue  of  the  Jacob ian  to  (7.3a)  are  A  =  -2000.5  and 
A„  =  -0.5  thus  yielding  a  stiffness  ratio  4001.   The  stiffness  matrices 
which  have  constant  entries  were  obtained  as 
-0.499875      0 

0      -0.499875 
2000.5    0 
i    0     2000.5 


(7.3b) 


n  - 


fi2  = 


/ 


Table  (7.3)  gives  the  numerical  results  for  the  explicit  scheme,  Lambert's 
scheme  and  the  blended  DIFSUB.   The  global  relative  error  for  the  blended 
DIFSUB  is  0.743400  D-5  while  the  global  relative  error  for  the  explicit 
scheme  is  0.5746777037  E-05.   The  computational  cost  and  function  evaluation 
on  IBM  360/75  for  the  blended  DIFSUB  are  respectively  $1.84  and  107  while 
that  of  the  explicit  scheme  is  $0.77  and  10,  respectively. 
The  global  errors  were  computed  as 


(7.3c) 


(7.3d) 


m 
e  =  max    E 
0<t<N  i=l 


yt  "  y(xt) 


1/2 


i  fl      i    i        i   , 

co  =  max   {1,   yQ,   y^...,  y^i 

l<i<m 
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Example  4  (see  Enright,  et  al.  [5] 

Van  der  Pol's  oscillator  (considered  in  Enright,  et  al.  [6]) 
(7.4)    y{  -  y2  y1(0)  =  2 

y»  =  5(l-y^)y2  -  Y;L    y2(0)  =  0 

0  £  X  £  1. 

Eigenvalues:  -0.067  and  -15  -+  5.7  and  -1.5  +  3.6  and  1.4  ->  2.4  ±  2.8i 
-+  -0.052  ±  8.8i  ■>  -2.0  -*•  9.5i  ■>  -5.9  ■>  4.5i  ■*  -2.0  and  -12  -v  0.050  and 
-15  +  1.1  and  -3.4. 

The  numerical  computation  was  effected  with  both  the  explicit 
and  the  implicit  schemes  using  uniform  mesh  sizes  h  =  0.2,  0.15,  0.1, 

0.05,  0.025  and  0.0125.   The  same  problem  was  solved  with  the  blended 

-3 
DIFSUB  using  an  initial  mesh  size  h  =  10   as  suggested  in  Enright, 

et  al.  [6]. 

From  Table  7.4  the  new  schemes  compare  favorably  with  the 

blended  DIFSUB. 
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TABLE  7.4. 

Numerical  Results  f 

or  Example  4. 

EXPLICIT 

No.  of 
Fn.  Eval. 

h 

yd) 

y(2) 

5 

0.2000 

1.8716065 

-0.14358810 

7 

0.1500 

1.8711045 

-0.14574713 

10 

0.1000 

1.8705973 

-0.14610294 

20 

0.0500 

1.8694380 

-0.14823599 

40 

0.0250 

1.8694389 

-0.14823587 

80 

0.0125 

1.8694388 

-0.14823588 

Blended  DIFSUB 

No.  of 
Fn.  Eval. 

h0=10"3 

e  =io-6 

y(D 
1.86944 

y(2) 
-0.148236 

IMPLICIT 

No.  of 
Fn.  Eval. 

h 

yd) 

y(2) 

- 

0.2000 

* 

* 

- 

0.1500 

* 

* 

42 

0.1000 

1.8693953 

-0.14824187 

81 

0.0500 

1.8694357 

-0.14823631 

161 

0.0250 

1.8694387 

-0.14823589 

321 

0.0125 

1.8694389 

-0.14823587 

*  No  convergence 
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Example  5 


(7.5a) 


/-icf5 


100 


\ 


0 


y'  = 


-5 


y(0)  = 


0  <  X  <  lOlT. 


1 


\-100   -10 

This  is  the  model  problem  (1.2)  with  £  =  10   and  w  =  100  whose 

theoretical  solution  is 

_in~5   sin  lOOx 
(7.5b)    y(x)  =  e  iU  X 

\  cos  lOOx 

The  numerical  results  in  Table  (7.5)  are  obtained  with  a  uniform 

mesh  size  h  =  tt/20  and  only  one  evaluation  of  the  stiffness/oscillatory 

parameters  was  obtained  and  given  by 


-10"5-1000i  o\ 


?'-10  5+1000i  o\ 


°1  =   1       -5       I'   "2  =  I       -5 

\0  -10  -10001/        \0     -10  +1000i 


TABLE  7.5.   Numerical  Results  for  Example  5. 
x  =  10  ,  uniform  h  =  tt/20 


X 

in12        Xt 
10        x     Tt+1 

12        2 
10±Z   x     T 

7T 

0.16918 

0.00949 

2tt 

0.02447 

0.02233 

3tt 

0.29342 

0.03504 

4tt 

0.61126 

0.04402 

5tt 

0.92906 

0.05704 

6tt 

1.2468 

0.07060 

7T7 

1.56451 

0.08478 

8tt 

0.09729 

0.09991 

9tt 

1.29056 

1.06429 

10tt 

1.60815 

0.12150 

t+1 


one  evaluation  of  oscillatory/stiffness  parameters 
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Example  6  (Stiefel  and  Bettis  [24]) 

We  finally  consider  the  nearly  periodic  initial  value  problem 
which  was  earlier  studied  by  Stiefel-Bettis  [24]  and  Lambert-Watson  [16]: 

(7.6a)    yM  +  y'  =  0.001eix,  y(l)  =  1,  y'(0)  =  0.99951,  y  e   R1 

whose  theoretical  solution  is 

y(x)  =  u(x)  +  iv(x) ,    u,  v  G  R 
(7.6b)    <(  u(x)  =  cos  x  +  0.0005x  sin  x 
j  v(x)  =  sin  x  -  0.0005  cos  x 

Equations  (7.6b)  represent  motion  on  a  perturbation  of  a  circular 

orbit  in  the  complex  plane  in  which  the  point  y(x)  spirals  slowly  outwards 

such  that  its  distance  from  the  origin  at  any  time  x  is  given  as 

(7.6c)    t(x)  =  /(u2(x)  +  v2(x)) 

The  initial  value  problem  (7.6a)  can  be  expressed  as 

|V  =  2y  Vo)  =  i 

V  =  -S  +  0.001  cos  x  ,   2y(0)  =  0 
(7.6d)    { 

|  V  =  y  ,   y(0)  =  o 

V  =  "3y  +  0.001  sin  x  ,   4y(0)  =  0.9995 

The  system  (7.6d)  was  solved  with  the  explicit  formulas  in  the 
range  0  <_  x  _<  40tt  which  corresponds  to  20  orbits  of  the  point  y(x).   The 
integration  was  performed  using  uniform  mesh  sized  h  =  tt/4,  tt/5,  tt/6,  tt/9 , 
and  tt/12.   Two  sets  of  numerical  results  were  generated — in  the  first  set, 
the  oscillatory/stiffness  parameters  are  obtained  once  at  the  first  step 
of  integration  while  in  the  second  set,  the  oscillatory  parameters  are 
evaluated  at  every  step  of  integration. 
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The  same  problem  was  solved  with  the  symmetric  multistep  method 
of  Lambert  and  Watson  [16]  as  well  as  the  StBrmer-Cowell  five-step 
multistep  formula  [24]  (both  of  order  6): 

(7.6e)    yt+5  -  2^  +  y^  -  JL  (18ft+5  +  209^  +  kt^   +  l«t+2 

,    -6ft+l+ft>- 

The  exact  distance  from  the  origin  x(x)  is  given  by  (7.6c) 

-jo    3  2 
and  the  approximate  distance  is  t  =  /(  y  +  y  )  at  x  =  40tt.   All  the 

solutions  generated  by  the  new  scheme  spiral  outward  in  agreement  with 

the  theoretical  solution  as  well  as  Lambert's  scheme  whilst  the  first 

three  values  generated  by  Stormer-Cowell  scheme  spiral  inward. 

From  the  Tables  (7.6a,  b  and  c)  which  respectively  show 

x,  |t(x)-t|  and  |y(x)-y|  =  /[  (  y(x)-  y)   +  (  y(x)-  y)  ] ,  we  see  that  despite 

the  fact  that  the  new  scheme  is  of  lower  order,  yet  it  is  more  accurate 

than  both  the  symmetric  multistep  method  as  well  as  the  five  order 

Stormer-Cowell  multistep  scheme. 
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TABLE   7.6a.      xf  =  40ir ,    t  (x  )    =  1.001972 


h 

Stbrmer-Cowell 

Symmetric 

Explicit 

one  evaluation 

repeated  eval. 

of  parameters 

of  parameters 

tt/4 

0.965645 

1.003067 

1.002311 

1.001972 

tt/5 

0.993734 

1.002217 

1.002205 

1.001972 

tt/6 

0.999596 

1.002047 

1.002140 

1.001972 

tt/9 

1.001829 

1.001978 

1.002050 

1.001972 

tt/12 

1.001953 

1.001973 

1.002016 

1.001972 

TABLE   7.6b, 


xf  =   40rr,   x  (xf )    =   1.001972 


h 

10  x|t(x)-t| 

109x|t(x)-t| 

Stormer-Cowell 

Symmetric 

Explicit 
one  evaluation  repeated  eval. 
of  parameters   of  parameters 

tt/4 

36  327 

1  095 

339 

204 

tt/5 

8  238 

245 

233 

66 

tt/6 

2  376 

75 

167 

26 

tt/9 

143 

6 

78 

3 

tt/12 

19 

1 

44 

0 

TABLE   7.6c.      xf  =   40tt,    u(x   )    =  1,   v(xf)    =   0.062832 


h 

106x  y(x J-yJ 

109x  y(xf)-y  | 

Stormer-Cowell 

Symmetric^ 

Expli* 
one  evaluation 
of  parameters 

:it 

repeated  eval. 
of  parameters 

TT/4 

48  014 

31  272 

389 

384 

tt/5 

13  136 

7  300 

252 

159 

tt/6 

4  494 

2  303 

176 

77 

tt/9 

405 

188 

79 

15 

tt/12 

73 

33 

45 

5 
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8.   CONCLUDING  REMARKS 

The  proposed  explicit  scheme  (2.3)  and  (2.9)  is  considered  to 
be  more  efficient  and  accurate  than  the  DIFSUB  and  the  blended  DIFSUB  for 
linear  stiff  systems  of  ordinary  differential  equations.   It  is  equally 
efficient  for  highly  oscillatory  systems  as  it  is  capable  of  admitting 
fairly  large  mesh  size  and  still  maintains  high  degree  of  accuracy.   The 
major  drawback  is  the  need  to  generate  higher  order  derivative  but 
automatic  generation  of  higher  order  derivatives  is  practicable  for  an 
extensive  range  of  problems. 
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